#!/usr/bin/env python
#coding=utf-8

import numpy as np
from pprint import pprint
from itertools import product
import math
from math import fabs
from numpy import array
import matplotlib.pyplot as plt
import os.path as path
from circleUtils import radius
import json
from circleUtils import edges


def getPointThroughput(x):
    """
    不能直接使用edges中的同名函数
    因为circleUtils默认的points与本地的points组织形式不同，前者最前面有四个物流园区点
    """
    assert type(x) is int
    return edges[x+4].sum() + edges.transpose()[x+4].sum()

points = []
def generateFirstLevelNodes(optimize = False):
    # 读入所有点
    global points
    with open("points.txt") as f:
        for line in f:
            x, y = line.split(',')
            points.append(np.array((float(x), float(y))))

    print "Points Number: ",len(points)

    # 找出平面上所有点两两之间距离小于R的集合
    candidates = candidatePointsPair(points)

    # 已知两个圆周上的点和半径可以确定一个圆，求出这些候选点集所对应的圆的圆心
    circles = []
    for p in candidates:
        circles.extend(getCenter(*p))

    # 将所有点本身也加入备选圆的集合中，以防止孤立点（无法在给定条件下与其它点构造出圆）出现
    circles.extend(points)

    print "Circle Number: ", len(circles)

    # 求出每一个圆所包含的点的集合
    index_points = zip(range(len(points)),points)


    print "Calculating point pairs..."
    circle_cover = []
    circle_cover_cache_file_name = ".circle_cover.cache" if not optimize else ".circle_cover.opt.cache"
    def load_circle_cover(circle_cover):
        if path.isfile(circle_cover_cache_file_name):
            print "load circle cover cache"
            with open(circle_cover_cache_file_name) as f:
                circle_cover.extend(json.load(f))
            return True
        return False

    if not load_circle_cover(circle_cover):
        for c in circles:
            circle_cover.append( (tuple(c), [x[0] for x in index_points if distance(c, x[1])<=radius+1e-3]) )
        with open(circle_cover_cache_file_name, "w") as f:
            json.dump(circle_cover, f)


    print "Candidate circles: ", len(circle_cover)

    # 按照圆覆盖的点的数量多少，每次选取一个能覆盖最多未被覆盖的点的圆，加入到所选择的集合当中。
    # 对于能添加的新点一样多的圆，应该选取本身包括的点更多的

    print "Calculating selected cover circles..."
    selected_circles = []

    selected_cache_file_name = ".selected.cache" if not optimize else ".selected.opt.cache"
    def load_selected_circle(selected_circles):
        if path.isfile(selected_cache_file_name):
            print "load selected circles cache"
            with open(selected_cache_file_name) as f:
                selected_circles.extend(json.load(f))
            return True
        return False


    coverred_points = set()

    if not load_selected_circle(selected_circles):
        while len(coverred_points) != len(points):
            c, s = sorted(circle_cover, key=lambda x:len( set(x[1]) -coverred_points ))[-1]
            coverred_points |= set(s)
            selected_circles.append(c)
        with open(selected_cache_file_name, "w") as f:
                json.dump(selected_circles, f)

    selected_circles = [tuple(x) for x in selected_circles] # 将每个元素转换成tuple以便于后续操作
    selected_circles_cover_dict = {tuple(x[0]):x[1] for x in circle_cover if tuple(x[0]) in selected_circles} # 以圆心为键，set为值

    # 优化：存在一些节点，去掉它们只会影响单独的一个区域，而这个区域本身能够通过一个更廉价的二级节点被包括到
    # 因此将其去掉并以一个二级节点代替
    if optimize:
        # selected_circles = [tuple(x) for x in selected_circles] # 将每个元素转换成tuple以便于后续操作
        # selected_circles_cover_dict = {tuple(x[0]):x[1] for x in circle_cover if tuple(x[0]) in selected_circles} # 以圆心为键，set为值

        whole_set = set(selected_circles)
        whole_points_set = set.union( *[set(selected_circles_cover_dict[x]) for x in selected_circles])

        selected_circles_set = set(selected_circles)


        for s in whole_set:
            remain = selected_circles_set - set([s])
            remain_sum_set = set.union(*[set(selected_circles_cover_dict[x]) for x in remain]) 
            remain_point_set = set.union(*[set(selected_circles_cover_dict[x]) for x in selected_circles_set]) - remain_sum_set
            if len(remain_point_set)==1:
                point = list(remain_point_set)[0]
                if getPointThroughput(point)>3000:
                    continue

                selected_circles_set.remove(s)
                
        selected_circles = list(selected_circles_set)
        selected_circles_cover_dict = {tuple(x[0]):x[1] for x in circle_cover if tuple(x[0]) in selected_circles} # 以圆心为键，set为值


    print "Selected circle number: ", len(selected_circles)

    return selected_circles_cover_dict

#暂时未使用该函数
def selectCircle():
        """
        选取一个能覆盖最多未被覆盖的点的圆
        对于能添加的新点一样多的圆，应该选取本身包括的点更多的
        """
        global circle_cover, coverred_points
        maximum_extra_points = 0
        total_points = 0
        result = None
        for c, s in circle_cover:
            extra = len(set(s) - coverred_points)
            if  extra > maximum_extra_points:
                maximum_extra_points = extra
                total_points = len(set(s))
                result = c,s

            if extra == maximum_extra_points and len(set(s))>total_points:
                total_points = len(set(s))
                result = c,s

        return result
def distance(a, b):
    return np.linalg.norm(a-b)


def candidatePointsPair(points):
    return [x for x in product(points, points) if not (x[0]==x[1]).all() and distance(*x)<2*radius]


def find_max(c, key):
    m = None
    result = None
    for x in c:
        if m is None:
            m = key(x)
            result = x
            continue
        if key(x) < m:
            m = key(x)
            result = x
            continue
    return result


def getCenter(p1, p2):
    """
    根据两个圆周上的点以及确定的半径，求圆心坐标
    """
    def sq(x):
        return float(x*x)
    a = p1[0]
    b = p1[1]
    x = p2[0]
    y = p2[1]
    R = S = radius
    d = math.sqrt(sq(math.fabs(a-x)) + sq(math.fabs(b-y)))
    assert fabs(d) > 1e-3
    assert d < (R+S) and d > (math.fabs(R-S))    

    A = (sq(R) - sq(S) + sq(d)) / (2 * d)
    h = math.sqrt(sq(R) - sq(A))

    x2 = x + A * (a-x)/d
    y2 = y + A * (b-y)/d

    x3 = x2 - h * ( b - y ) / d
    y3 = y2 + h * ( a - x ) / d

    x4 = x2 + h * (b - y) / d
    y4 = y2 - h * (a - x) / d

    return ((x3,y3),(x4,y4))





if __name__ == "__main__":
    
    selected_circles_dict = generateFirstLevelNodes(True)
    # pprint(selected_circles_dict)
    selected_circles = selected_circles_dict.keys()
    def drawCircle(c):
        a=c[0]
        b=c[1]
        r= radius
        theta = np.arange(0, 2*np.pi, 0.01)
        x = a + r * np.cos(theta)
        y = b + r * np.sin(theta)
        plt.plot(x,y)

    fig = plt.figure() 
    fig.set_size_inches(18.5, 10.5)


    # big_points = []
    # for k, p in enumerate(points):
    #     if edges[k].sum() + edges.transpose()[k].sum()>3000:
    #         big_points.append(p)
    
    x = [t1[0] for t1 in points]
    y = [t1[1] for t1 in points]
    plt.scatter(x,y)

    # x = [t1[0] for t1 in big_points]
    # y = [t1[1] for t1 in big_points]
    # plt.scatter(x,y,marker='x')

    for c in selected_circles:
        drawCircle(c)

    plt.show()

